【直播】我的基因组72：把基因检测芯片数据转为vcf格式

原创 2017-04-27 jimmy 生信技能树

这个需求比较少见，主要是因为有很多朋友都做了基因检测芯片数据，而芯片检测的结果只有4列，分别是dbSNP数据库ID号，染色体，坐标，还有基因型。如下：

但是呢，大部分的基因检测结果注释都是基于vcf文件的，vcf文件的详细介绍，我们以前讲过，就是

这10列。

要想把基因检测芯片数据转为vcf格式就需要在充分理解vcf的基础上面再增加几个信息。

因为基因芯片的结果里面没有参考碱基是什么的信息，只有基因型，所以我们没办法判断纯合杂合或者突变。

至于 QUAL，这里统一给100，filter统一给一个占位符即可， INFO就给一个测序深度， FORMAT就给 GT:DP:RO:QR:AO信息吧，具体介绍如下：

##INFO=<ID=DP,Number=1,Type=Integer,Description="Total read depth at the locus">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##FORMAT=<ID=RO,Number=1,Type=Integer,Description="Reference allele observation count">
##FORMAT=<ID=AO,Number=A,Type=Integer,Description="Alternate allele observation count">

这里我们还是借用dbSNP数据库文件，如下：

这个文件我以前讲过：

那么很简单的一个perl程序就可以达成这个转换啦：

open FH,"wegene.txt";
while(<FH>){
chomp;
@F=split;
next if /^#/;
$h{$F[0]}=$F[3] if /[ATCG]{2}/;
}
close FH;
open FH,"/home/jianmingzeng/annotation/variation/human/dbSNP/dbsnp.pos";
open OUT,">wegene.vcf";
print OUT '##INFO=<ID=DP,Number=1,Type=Integer,Description="Total read depth at the locus">'."\n";
print OUT '##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">'."\n";
print OUT '##FORMAT=<ID=RO,Number=1,Type=Integer,Description="Reference allele observation count">'."\n";
print OUT '##FORMAT=<ID=AO,Number=A,Type=Integer,Description="Alternate allele observation count">'."\n";
print OUT '#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT unknown'."\n";
while(<FH>){
chomp;
@F=split;
next unless exists $h{$F[2]};
next if $h{$F[2]} eq $F[3].$F[3];
$gt=$h{$F[2]} eq substr($h{$F[2]},0,1)x2 ?"1/1":"0/1";
$alt=substr($h{$F[2]},0,1) eq $F[3] ? substr($h{$F[2]},1,1): substr($h{$F[2]},0,1);
$ro_po=$gt eq '1/1' ? "0:100" : "50:50" ;
#print "$_\t$gt\t$h{$F[2]}\n";
print OUT "$F[0]\t$F[1]\t$F[2]\t$F[3]\t$alt\t100\t.\tDP=100\tGT:DP:RO:AO\t$gt:100:$ro_po\n";
}
close FH;
close OUT;

运行完毕就可以打开我们转换好的vcf文件，如下所示：